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Abstract 

We examine axino dark matter in the regime of a low reheating temperature, Tr, after inflation and 
taking into account that reheating is a non-instantaneous process. This can have a significant effect on the 
dark matter abundance, mainly due to entropy production in inflaton decays. We study both thermal and 
non-thermal production of axinos in the framework of the MSSM with ten free parameters. We identify the 
ranges of the axino mass and the reheating temperature allowed by the LHC and other particle physics data 
in different models of axino interactions. We confront these limits with cosmological constraints coming the 
observed dark matter density, large structures formation and big bang nucleosynthesis. We find a number of 
differences in the phenomenologically acceptable values of the axino mass mg and the reheating temperature 
relative to previous studies. In particular, an upper bound on mg becomes dependent on Tr, reaching a 
maximum value at Tr ~ 10 2 GeV. If the lightest ordinary supersymmetric particle is a wino or a higgsino, 
we obtain a lower limit of approximately 10 GeV for the reheating temperature. We demonstrate also that 
entropy production during reheating affects the maximum allowed axino mass and lowest values of the 
reheating temperature. 
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1 Introduction 


A variety of cosmological data points to the existence of a non-luminous component of the matter in the 
Universe, referred to as dark matter (DM). In spite of decades-long efforts aiming at detecting DM particles 
directly, they have so far remained elusive. Indirect probes, such as the impact of DM on the formation and 
evolution of large structures in the Universe or on the spectrum of the temperature fluctuations in the cosmic 
microwave background provide powerful tools to study DM. However, such probes involve only some properties 
of DM particles, such as their (time-dependent) energy density and free-streaming length, hence the masses of 
proposed DM candidates span thirty orders of magnitude while their cross-sections for scattering with known 
particles - forty orders of magnitude (see, e.g., 1) for a recent review). 

There are many examples of DM candidates which were not introduced ad hoc , but possess a sound theoretical 
motivation. They include bosonic fields in coherent motion, such as axions, weakly interacting massive particles 
(WIMPs), such as the lightest neutralino of the Minimal Supersymmetric Standard Model (MSSM), or extremely 
weakly interacting particles (EWIMPs), such as gravitinos and axinos. 

The axion was originally introduced as a solution to the so-called strong CP problem. The smallness of the 
electric dipole moment of the neutron can be understood if there is a C/(1)pq symmetry, referred to as Peccei- 
Quinn (PQ) symmetry [2113], which is spontaneously broken at a scale f a . The very light pseudo-Goldstone 
boson a associated with this symmetry, called an axion, couples to the gluon anomaly Hj (see also, e.g., [5] for 
a review), 
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where a s = g 2 /47 t and g s is the strong coupling constant. Various cosmological considerations suggest that f a 
lies between approximately 10 9 GeV and 10 12 GeV )6j. 

Two broad frameworks for field theory models of axions have been considered in the literature. In the 
Kim-Shifman-Vainshtein-Zakharov (KSVZ) model [7, 8], the gluon anomaly term 0 is induced through a 
loop of very heavy singlet quarks carrying PQ charges, while in the Dine-Fischler-Srednicki-Zhitnitsky (DFSZ) 
model E1HD3 the PQ charges are assigned to the Standard Model fields. These assumptions are sufficient to 
determine unambiguously axion interactions with the Standard Model fields that are relevant for cosmology, 
though specific implementations of these models can be rather complicated SHU- The Lagrangian describing 
axion interactions with Standard Model particles is given in Appendix A. 

In supersymmetric models [ T2] the axion resides in a chiral multiplet with a fermionic superpartner - 
the axino a - whose interactions with the Standard Model particles are related to those of the axion through 
supersymmetry. See BUS] for a review, relevant formulae are shown in Appendix A. The properties of the axino 
in models with softly broken supersymmetry can be relevant for cosmology [Musm especially when it is the 
lightest supersymmetric particle (LSP) and constitutes cold DM I?, IB]- Particular scenarios of supersymmetry 
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breaking give predictions for the axino mass HS1 Eol El]; however, due to a strong model dependence and the 
absence of an universally accepted scheme of supersymmetry breaking, we adopt a phenomenological approach 
and treat the axino mass as a free parameter. 

The mechanisms for axino generation in the post-inflationary Universe (see |22l for a review) are analogous 
to those for the gravitino: thermal production (TP) from scatterings and decays of other particles in thermal 
equilibrium and non-thermal production (NTP) from out-of-equilibrium decays of heavier particles. Although 
detailed predictions for the axino abundance from TP, Y" 5 TP , strongly depend on the model of axion interactions, 
it is inversely proportional to f 2 and does not depend on the axino mass. In this respect, the axino differs 
significantly from the gravitino: as the latter contains a goldstino component related to broken supersymmetry, 
the gravitino abundance from TP is inversely proportional to Mp and to the square of the gravitino mass. A 
notable difference between the axino interaction models is that in KSVZ models is proportional to the 
reheating temperature Tp defined in terms of the inflaton decay rate = Tf\/g*{Tp)/90 (T^/Mpi), while in 
DFSZ models U 5 TP does not depend on Tp. This can be attributed to the fact that in DFSZ models axinos 
are mainly produced from decays of thermal particles, while in KSVZ models the main source of axinos are 
scatterings of strongly interacting particles. Existing numerical analyses [23JEU suggest that, in order not to 
overdose the Universe in the KSVZ scheme, Tp should be at most a few orders of magnitude larger than the 
electroweak scale. 

The contribution of the NTP to the axino abundance comes from decays of the lightest ordinary super- 
symmetric particles (LOSP) which underwent freeze-out, hence V~ NTP = Vlospi or in terms of cosmological 
parameters, 
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Depending on the masses of the axino and the LOSP, the LOSP lifetime may be so long (above seconds) that the 
highly energetic LOSP decay products can affect successful predictions of the big bang nucleosynthesis (BBN). 

Although there has been remarkable progress in the understanding of axino physics and cosmology in the 
last two decades, the latest LHC data and the increasing precision of reconstructing the history of the Universe 
motivate extending the existing analyses in two ways. Previous studies of axino DM focused on axino interac¬ 
tions, while some specific assumption about the MSSM mass spectra were made for simplicity. Also, the low Tp 
regime was studied assuming instantaneous reheating and some fixed typical values of the abundance of axino 
DM originating from NTP. Given many constraints that the LHC data put on the MSSM, it is worthwhile to 
ask what are the allowed and excluded ranges of the parameters of the general MSSM with axino DM. It is also 
known that the altered expansion rate of the Universe during reheating may result in a significant change of 

DM abundance E3E1EI3EBI2H1I1BI3III32]- 

In this paper, we address the two issues mentioned above. We identify the phenomenologically viable pa¬ 
rameter ranges of the MSSM with axino DM and calculate accurately the axino abundance [33] not only during 
the radiation dominated (RD) period, but also during the phase of reheating after inflation [31] , The paper is 
organized as follows. In Section[2j we examine how non-instantaneous reheating affects the predictions for axino 
DM. In Section [3] we show the results of our numerical study of a 10-parameter version of phenomenological 
MSSM (plOMSSM) and identify the ranges of values of the axino DM mass and the reheating temperature that 
are consistent with data including large scale structure formation and big bang nucleosynthesis (BBN) con¬ 
straints. We present the conclusions in Section [4] A number of technical issues are relegated to the Appendices. 
Appendix A summarizes the interaction Lagrangians of the Standard Model or MSSM fields with the axion and 
the axino. Appendix B presents some details of our calculation of TP of axinos and Appendix C describes some 
phase space integrals necessary to address the free-streaming length of axinos from NTP. In Appendix D, we 
describe the details of our numerical procedure. 
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Figure 1: Left panel: The ratio of the axino abundances in the KSVZ model obtained with the assumption 
of non-instantaneous (yi rp ’ non-lnst ) anc j instantaneous (y- TP,lnst ) reheating as a function of Trd- The gluino 
and squark masses are set to 1 TeV. Right panel: Predictions for TP of axinos resulting from scatterings and 
squark decays for 10 2 GeV < Trd < 10 6 GeV; results obtained for instantaneous reheating are also shown. 

2 Axino with low reheating temperature — general remarks 

2.1 Thermal production of axinos 

Our first aim is to understand the effects of low Tr and of non-instantaneous reheating on axino abundance. 
Therefore, we will first discuss TP of axinos, focusing mainly on the KSVZ model. For concreteness, we will 
present the results obtained for the gluino and squark masses of rrig = rriq = 1 TeV, and we shall assume that 
the C a ww = 0 (see Appendix A for details regarding parameters and notations), unless indicated otherwise. 
We shall relax this assumption about the MSSM mass spectrum later. 

In the high-Tfj regime, the scatterings associated with the SU( 3) c group dominate axino TP0 Were the 
reheating process instantaneous, we could identify the reheating temperature Tr with the temperature Trd at 
which the radiation dominated epoch began. However, with the RD epoch having been preceded by a reheating 
epoch, during which the energy density of the oscillating and decaying inflaton dominated the Universe, there 
is an additional contribution to Y/f p originating from a modified relation a(t) ~ T -8 / 3 * * between the scale factor 
of the Universe a(t) and the temperature T. A straightforward, but technical and rather involved calculation 
presented in Appendix B leads to the conclusion that this additional contribution is about 1/6 of the standard 
high Tr result. Loosely speaking, the existence of a reheating phase before the RD epoch effectively extends 
the period of relic production and available temperature range - and therefore the axino abundance increases. 

One comment is in order here. In Ref. |3()j , which studied the same situation as described above, a constant 
reduction of YA p by a factor of about 1/4 was reported. This apparent difference results from a different choice 
of a reference point. In Ref. m the results obtained for instantaneous and non-instantaneous reheating were 
compared for the same value of Tr, while we believe that it is more appropriate to compare both abundances 
at the same value of Tr Dj which (in contrast to Tr being a convenient shorthand notation for the inflaton 
decay rate) has a clear physical meaning of the temperature that marks the beginning of the standard radiation 
dominated epoch. As shown in Appendix B, we find an approximate relation Trd ~ Tr/ 2, so we can use these 

2 In literature, there exist different prescriptions for treating the infrared divergence in relevant scattering cross-sections (see, 

e.g., [T] . We conclude that there currently remains a factor of a few uncertainty in the thermal yield of axinos at high Tr. In 

our numerical analysis, we will use the effective mass approximation m ■ However, as we will focus on the case of low Tr, this 

uncertainty will be of no consequence for our main results; see Section [3| 
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Figure 2: Predictions for TP of axinos for different choices of the MSSM mass spectra, as indicated in the plots. 
Non-colour contribution is negligible here. 


two temperatures interchangeably when referring to orders of magnitude. 

For intermediate values 10 2 GeV ;$ Trd ~ 10 4 GeV, there is a phase space suppression of the scattering 
terms associated with TeV-scale superparticles. For a given value of Trd this effect is smaller in the case 
of non-instantaneous reheating, as the additional contribution to Y? p stabilizes the result and the ratio of 
the abundances calculated for non-instantaneous and instantaneous reheating becomes slightly larger than 
7/6 w 1.17. This can be seen in the left panel of Fig. [T] 

For Trd smaller than the masses of strongly interacting particles (herein 1 TeV) the ratio falls, as the 
contribution to axino TP from scatterings becomes small with respect to the contribution from squark and 
gluino decays. This is because, unlike scatterings, the decays depend rather weakly on the details of reheating 
as shown in the right panel of Fig. [T] In principle, larger temperature values attainable with non-instantaneous 
reheating may lead to a larger equilibrium number density of decaying squarks or gluinos and therefore to an 
increase of F~ rp , but this effect is less important than the phase space suppression in the case of scatterings. 
Hence, as long as Y a TP is dominated by decays, the ratio in the left panel of Fig. [l] falls below 1.17. Fig.[2]shows 
examples of the dependence of axino TP on the gluino and the squark masses. 

The increase of the ratio of the abundances for the instantaneous and non-instantaneous reheating scenarios 
can be sizable only for Xrd < 100 GeV, if the parameter C a yy parametrizing the coupling between the axino 
and the U(l)y gauge bosons and gauginos is equal to zero or if bino-like neutralino is heavier than about 
500 GeV. Otherwise, axino TP is at low temperatures dominated by the decays of the light bino which is again 
rather insensitive to the details of reheating if Trd is comparable to the bino mass. 

However, for such low values of Trd, even for vanishing C a yy and/or heavy bino, TP often gives an 
abundance which is a very small fraction of that required for DM, so details of reheating are irrelevant in that 
case. One should also remember that for very low values of Trd axinos are decoupled from kinetic but not from 
chemical equilibrium and their distribution function can differ from that for the equilibrium case [34] . but for 
axinos this only happens for the DM abundance dominated by NTP, so this effect has no consequences for our 
study. 

We can see that the main effect of non-instantaneous reheating on the axino TP originates from a modified 
contribution from scattering. Therefore, the details of reheating play a minor role in DFSZ models, as in this 
case TP is typically dominated by thermal higgsino decays [231 i25 • 
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Figure 3: Cosmological parameters D 5 / 1 2 resulting from TP and NTP of KSVZ axinos for different values of 
the reheating temperature Tr. Gluino and squark masses are set to 1 TeV, while the bino masses was set to 
100 GeV with the rest of the MSSM spectrum arranged so that the high-7^ bino abundance is il^h 2 = 10 3 4 
(left panel) and to 1 TeV with the rest of the MSSM spectrum arranged so that the high-T# bino abundance is 
flgh 2 = 0.1 (right panel). 

2.2 Non-thermal production of axinos 

For low enough values of T R d axino TP is highly suppressed and it is the non-thermal contribution to the axino 
relic density that dominates. However, NTP is also affected by an additional entropy production during the 
reheating period provided that T RD is low enough so that the LOSPs freeze out before the RD epoch (see, 
e.g., a recent discussion about similar issue in the case of gravitino DM [3TJ)|^] This results in a suppression 
of Dlosp^ 2 below the value flLOSP^ 2 (high T R d) obtained in the standard cosmological scenario where the 
LOSP freeze-out occurs in the RD epoch. Consequently, for a fixed ma/mLOSP ratio in Eq. ([ 2 ]), Qf TP h 2 also 
decreases. The lower the reheating temperature, the longer the period between the LOSP freeze-out and the 
beginning of the RD epoch, so the LOSPs are effectively diluted. As a result r2j TP /i 2 decreases with T RD . This 
is illustrated in Fig. [ 3 ] where we show both TP and NTP contributions to for two selected SUSY spectra 
with bino LOSP mass equal to 100 GeV and ~ 1 TeV, while squark and gluino masses are set to 1 TeV. The 
rest of the SUSY spectrum is in both cases chosen such as to obtain LOSP yield at freeze-out corresponding to 
Og/i 2 (high T r ) = 0.1 or 10 4 . 

2.3 BBN and WDM constraints 

The BBN constraints for the axino can be analyzed similarly to those for the gravitino (see e.g. [33 EE! 022 
HOI HU HU mi HU); they are typically mild. This is because the lifetime of the LOSP decaying to the axino 
usually hardly exceeds 0.1 sec., unless one considers very light LOSP, allows for a very strong mass degeneracy, 
rria ~ tolosp or considers a LOSP whose 2-body decays to axinos are forbidden (i.e. bino with C a YY — 0). 
However, if the neutralinos are too light, they have a very small relic abundance, which is further suppressed by 

3 The most conservative lower bounds on the reheating temperature from BBN can be derived in terms of a successful neutrino 
thermalization for electromagnetic energy emission, T R > 0.7 MeV and T R > 4 — 5MeV for weak-scale parent particles in terms of 

n — p conversion due to hadron emissions (35i 361 . 
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a small value of Tr, and there are simply too few decaying LOSPs around to put primordial nucleosynthesis in 
danger (in this case, the correct axino DM abundance must originate from TP). The only exception is a scenario 
with a very light axino and a very light bino LOSP with extremely small annihilation rate, suppressed by large 
slepton masses. In this case, the relic abundance of bino LOSPs can be very large and the correct axino DM 
relic density is obtained from NTP by a suppression with a small axino mass. 

A more dramatic change appears when we set C a YY = 0, thereby assuming that the axino interacts only 
with the SU(3) C gauge sector. For bino LOSP, the 2-body bino decay into axino and a photon or Z boson 
are then disallowed and the dominant bino decay channel B —► qqa involves an effective axino-quark-squark 
interaction vertex discussed in |45j . Applying the formulae obtained therein, we can estimate the bino lifetime 
as (see also Appendix C) 



With a longer LOSP lifetime, the BBN constraints are significantly more severe. 

For small enough values of mg, axinos may become warm dark matter (WDM) with a free streaming length 
which is too large to account for observed structures in the Universe]^] Thermally and non-thermally produced 
axinos have different rms velocities and hence these constraints from structure formation depends on the fraction 
of the WDM axinos in DM (for a discussion about non-thermally produced WDM see, e.g., g?] 05]). We use 
the 95% CL exclusion limits from Ref. [49]. For TP domination these limits translate to mg < 5keV, while 
the case of NTP domination and the mixed case require a more careful analysis In particular, for C u yy = 0 the 
calculation of the rms axino velocity requires an integration over 3- and 4-particle phase space (bino and stau 
LOSPs, respectively); we provide technical details of this computation in Appendix C. 


3 Axino with low reheating temperature in the MSSM 

Having discussed the predictions for TP and NTP of axinos for different values of the reheating temperature 
T r and for different patterns of the MSSM mass spectra, we are now ready to present and discuss the results 
of our extensive numerical scan. There are shown in Fig. [I] for f a = 5 x 10 9 GeV and 10 11 GeV. In comparison 
with the results of previous studies, we find a number of differences in the shape of the allowed region of the 
axino mass and the reheating temperature. We discuss them separately for the bino LOSP and for the wino 
and higgsino LOSP. 

3.1 Bino LOSP 

Firstly, similarly to Ref. [33] we find an upper bound on Tr and a lower bound on mj, beyond which axino DM 
becomes too warm. For Tr < 10 6 GeV the region allowed by the BBN or WDM constraints extends to smaller 
values of mg ~ 10~ 4 GeV than in previous analyses. For Tr < 10 3 GeV these points correspond to dominant 
NTP from a very large relic density of bino LOSP with a very small annihilation cross-section dominated by a 
t-channel exchange of multi-TeV sleptons. Such points, are however excluded by both BBN constraints (bino 
LOSP has a long lifetime and a sizable hadronic branching fraction, cf. [5D]) and by WDM constraints (with 
dominant NTP, one needs mg > 30MeV). Eventually, the lower bounds for mg at fixed Tr are similar to those 
obtained in [35], but in our analysis they result from incorporating additional cosmological constraints. 

Secondly, looking at the rightmost part of each plot in Fig. [4] we see an upper bound on mg. In this 
respect our result visually resembles that of Ref. [33], but this similarity follows from very different physical 
assumptions. In Ref. [33] three typical choices of fixed axino abundance from NTP, Y a NTP , were considered, 
which led to upper bounds on axino mass for Vg NTP > 0. When NTP was the dominant source of axino DM, 
these bounds on mg did not depend on Tr. In contrast, in our analysis we obviously do not make such a 
restrictive assumption and the maximum value of mg that we obtain is related to the maximum value of the 

4 For alternative, cosmologically viable scenario with axino WDM from late-time saxion decays see ED- 
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Figure 4: Results of the numerical scan for bino LOSP projected onto the (to^, Tr) plane for the KSVZ axino 
with f a = 5 x 10 9 (left panel) and f a = 1 x 10 11 (right panel). The area bounded by thick solid (partially 
shaded) lines denotes acceptable regions for which MSSM parameters consistent with all the bounds listed in 
Appendix D, including DM energy density, could be found. Inside these regions, we marked in red the regions 
excluded by axino dark matter being too warm (WDM constraints). Dashed vertical lines denote the lower 
bounds on the axino mass coming from WDM constraints for TP only. Regions excluded by BBN constraints 
are either dashed ( C a YY =8/3) or marked with a dash-dotted lines (C a YY= 0). 


LOSP mass. More precisely, the largest value of ma is obtained for Tr ~ 0(lO 2 GeV), corresponding to the 
freeze-out temperature of the bino LOSP. This maximum ma corresponds to the largest available values of the 
bino LOSP mass and the largest available slepton masses (cf. Table |T]) . These two features have the same effect: 
one expects a larger relic density for heavy particles and there is also a big suppression of the annihilation 
cross-section due to large masses of the intermediate particles. We show these aspects of our results in the left 
panel of Fig. [5j where we show how the allowed region changes with different assumptions about the largest 
possible bino and slepton masses. 

Furthermore, we see a decrease of the maximum allowed iria for Tr falling below 0(lO 2 GeV). This can 
be understood taking into account that during reheating there is entropy production due to inflaton decays. 
As a consequence, the LOSP relic density becomes suppressed, if Tr falls significantly below the freeze-out 
temperature. For a given value of Tr the suppression is stronger for binos with larger masses. This confines 
the allowed region to smaller values of bino mass and, as ma < mg, to smaller values of axino mass. 

We also note that, for Tr < 10 4 GeV the calculations for TP of axinos cannot be fully trusted, as the SU(3) C 
gauge coupling becomes large (see T3] and references therein). This poses no problem for Tr < 10 2 GeV, as 
NTP of axinos is then dominant, but in the window 10 2 GeV < Tr < 10 4 GeV the upper bounds on rria or Tr 
should be treated as approximate. Changing fa mostly leads to a shift of the allowed region in the (to 5 ,Tr) 
plane, as shown in the right panel of Fig. [4] 

The difference between regions of the (mj, Tr) plane allowed in the KSVZ and DFSZ models is presented in 
the right panel of Fig. [ 5 J For large values of Tr > 10 5 GeV it can be traced back to the fact that ttj p h 2 scales 
as TOsTr and mg for the KSVZ and DSFZ models, respectively [231124j . For smaller values of Tr the bulk 
of the allowed region is very similar for both models. For a fixed ma ~ 1 MeV, there is still a small difference 
between the largest allowed Tr corresponding to the TP dominance. This results from different sources of 
thermal contributions to axino DM from decays. In KSVZ models the most important decays are those of 
colored particles; these are loop-suppressed with respect to decays of neutralinos with a non-negligible higgsino 
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Figure 5: The left panel is a blow-up of the bottom-right corner of the right panel of Fig. [4j boundaries of the 
allowed region corresponding to different ranges of the bino mass, scanned up to 1, 3 and 5 TeV, and different 
ranges of the stau mass, scanned up to 5 and 15 TeV are shown (cf. Table [l] for mass range dependencies). 
The right panel shows the results for the DFSZ model with f a = lx 10 11 GeV; for comparison the contours 
corresponding to the allowed region in the KSVZ model (right panel of Fig. [4]) are also shown. 


fraction in DFSZ models. Therefore, in DFSZ models TP of axinos from decays is much more efficient which 
may lead to overproduction of axino DM for values of ma and Xr corresponding to the correct axino DM density 
within the KSVZ scheme. 

3.2 Wino, higgsino and stau LOSP 

In Fig. [b] we show the allowed regions of the (ma,!#) plane for the wino and higgsino LOSP and in Fig. [7] the 
same for the stau LOSP, with constraints imposed in the same way as described in Section |3.1[ We again show 
the results for two values of f a = 5 x 10 9 GeV and 10 11 GeV. We restrict our analysis to the KSVZ model, since 
we have seen that for low Tr the predictions of the KSVZ an DFSZ models do not differ very much. 

The main difference with respect to the bino LOSP case (Fig. [4]) is that Tr is now bounded from below. 
This is because for winos, higgsinos and staus (unlike for binos) the masses of the states mediating annihilations 
cannot be much larger than the masses of annihilating particles. Hence, the annihilation cross-section of 
these particles cannot be made very small by increasing the mass of the intermediate states. Nonetheless, by 
varying MSSM parameters, we find examples of axino DM and the higgsino or wino or stau LOSP in which the 
observed DM abundance originates mainly from TP (upper parts of the allowed regions) or NTP (lower parts 
of the allowed regions). This may seem at odds with a recent analysis of axino DM in the DFSZ model with the 
higgsino LOSP [32] which studied both TP (freeze-in of axinos) and NTP (LOSP freeze-out), and concluded 
that TP always dominates. This apparent difference can be traced to the fact that in Ref. [32] two examples 
of the MSSM mass spectrum are studied, while our numerical analysis extends to larger values of the masses of 
the supersymmetric particles and thus allows larger LOSP relic abundances. This in turn can give rise to the 
observed axino DM abundance via NTP even with additional entropy production from inflaton decays. 

In the case of the stau LOSP one can obtain significant contribution to the axino relic density from TP 
also for low values of Tr ~ 10 GeV. It is due to decays of light stau LOSPs being still in thermal equilibrium. 
The mass of the axino is then typically significantly lower than Wf 1 , which suppresses the NTP contribution. 
However, one needs to take into account that such a scenario is constrained by the LHC searches for direct 











Figure 6 : The same as in 


Fig. 0 but for higgsino LOSP (left panel) and for wino LOSP (right panel). 


production of staus with missing energy m- When treating this we calculate the relevant production cross 
section with MadGraph5_aMC@NL0 [52| . 

3.3 Direct and cascade decays of the inflaton field 

In our analysis so far, we have made an assumption that there are no direct and cascade decays of the inflaton 
field to axino DM particles. Here, we would like to study the impact of such decays on the allowed ranges of 
axino mass and reheating temperature, following the model-independent approach used in End]. 

Our results are presented in Fig. [7] and [ 8 ] where we plot the allowed regions for the stau, bino and higgsino 
LOSP for different values of the dimensionless parameter 77 = b ■ (100 TeV/m^), where b is an average number 
of axino DM particles per inflaton decay and is the inflaton mass at the minimum of the potential. As 
inflaton decays provide an additional non-thermal component of axino DM (see a recent discussion in the case 
of gravitino m), the allowed region becomes extended towards smaller values of Tr at largest allowed values 
of m„. This is because the additional NTP from inflaton decays allows for a smaller contribution to axino DM 
density originating from LOSP decays, hence - for a fixed set of the MSSM parameters - for a smaller Tr and 
a larger suppression of LOSP abundance by entropy production in inflaton decays. 


4 Conclusions and outlook 

In this paper we have studied the impact of a low reheating temperature Tr on thermal and non-thermal 
production of axino DM, taking into account the non-instantaneous nature of the reheating process. We also 
extended previous studies by analyzing wide ranges of phenomenologically acceptable parameters of the 10 - 
parameter version of phenomenological MSSM instead of presenting the results for a single typical parameter 
choice. Comparing our results with previous works, we found a number of differences in the allowed ranges of 
the axino mass rria and the reheating temperature. In particular, depending on the choice of the axion model 
and the choice of the MSSM parameters, we showed that BBN constraints can exclude large portions of the 
parameter space corresponding mainly to non-thermal production of axino DM relevant for low Tr. We also 
demonstrated how entropy production during reheating affects the upper limits on the axino mass for a given 
range of the MSSM parameters. 
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Figure 7: The same as in Fig. [4j but for stau LOSP. Additionally, regions allowed with direct and cascade 
decays of the inflaton field (see Section 3.3) are shown. 




Figure 8: The impact of direct and cascade decays of the inflaton on the allowed regions of (ms, Tr) plane for bino 
LOSP (left panel) and higgsino LOSP (right panel) in terms of the dimensionless parameter ?/ = b -(100 TeV/m^) 
defined in the text. 
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The are a few directions in which the analysis presented herein could be extended. In our work, we relaxed 
the simplified assumption of instantaneous reheating, we still required that the maximum temperature during 
the inflaton-dominated period was sufficiently larger than the reheating temperature that the supersymmetric 
particles could reach thermal equilibrium. Although it is a realistic requirement, one can also envision scenarios 
with a very small energy density at the end of inflation. Additionally, it has recently been noted that the 
maximum temperature during reheating may not be as large as previously estimated [53]. Such cases are not 
included in our study and we leave them for future work. 

An issue which requires some care in models with low Tr is the origin of the primordial baryon asymmetry. 
Since we work in a supersymmetric setup, the Affleck-Dine mechanism [541 (see, e.g., [ 55] for a review) is a 
feasible options, though a detailed discussion is beyond the scope of our study. 
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Appendix A: Axion and axino interactions 


The effective axion interaction Lagrangian after integrating out all heavy PQ charged fields can be written, to 
the lowest order terms in 1 //„, as 


/^eff 
^a, int 


Cl ^—7—^ 75 g — Yi q (q L in qu e lC2a/fa + h.c.) 

Ja 


+ 


C3 


327T 2 f c 


aGG- 


C a WW 
327T 2 f a 


aWW 


CgYY 
327T 2 f a 


a B B + Uiept O 


(4) 


where (following a partial integration over on-shell quark fields) the ci term can be reabsorbed into the ci term. 
The KSVZ case can be identified with ci = 0, C 2 = 0, C 3 7 ^ 0, while the DFSZ one with ci = 0, C 3 = 0, C 2 7 ^ 0. 
General axion models can have both C 2 7 ^ 0 and C 3 7 ^ 0. C a ww and C a yY are model-dependent parameters 
that correspond to axino-gaugino-gauge boson anomaly interactions for the U(l)y and the SU( 2)^ groups, 
respectively. 

In a supersymmetrized version of an axion model m\mm the real scalar axion field a resides in a chiral 
supermultiplet since it is a gauge singlet. The other members of the axion supermultiplet are the fermionic 
superpartner axino a and the real scalar field sa xion s that provides a remaining bosonic degree of freedom 
on-shell. 

The interaction Lagrangian for the axion supermultiplet can be obtained by supersymmetrizing Eq. Q. In 
particular, the axino-gaugino-gauge boson and the axino-gaugino-sfermion-sfermion interaction terms are given 

by nans] 


/*eff 


175 [lb 71 9 b G% + £j~J~9 a ^9s r T a q 


167T f a 

5 75 [r ’ 71 + £r a ~ dWa E/ ° 52 f ~*° T “ ~ fD 

+ j Q i 6?/7 5 75 [r ’ 71 e / 5y f ~* Qy 


(5) 


where fn and / denote sfermions carrying non-zero T 3 and Y. respectively. 

A generic form of interactions between the axion and matter supermultiplets was considered in [56]. In 
particular, it was pointed out that, for upq > T > M$, where is the mass of the heaviest PQ-charged and 
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gauge-charged supermultiplet $, the axino-gaugino-gauge boson interaction term is suppressed by M|/T 2 . This 
is particularly important for the DFSZ axino, where <f> corresponds to the Higgs supermultiplets and therefore 
M$ = p (the higgsino mass). The dominant contribution to axino TP is then associated with a higgsino decay 
to the axino and the Higgs boson that is described by [551 Ell E3] 

-^ajrFSZ 3 c h — o, [Hd H u + H u Hj] + h.c. (6) 

Ja 

Appendix B: Calculation of axino TP with low reheating temperature 

In scenarios with non-instantaneous reheating in calculating U 5 TP one has to take into account a modified 
expansion rate of the Universe. Below we briefly describe the methodology that can be used to calculate the 
axino TP yield. 


Axino TP yield with non-instantaneous reheating. It results in a modification of temperature depen¬ 
dence on the scale factor T(a). The Boltzmann equation can then be written as 

~dT da = H ( Sscat + Sdec ) ’ (7) 

where X a = a 3 ns.. In that case the present-day axino abundance can be written as: 


Ha, 0 = 


1 


S(Aq Jt 0 


IX- t 


dlnT\-i A 2 


dA ) 


II (^scat 


Udec 


( 8 ) 


where T up corresponds to an effective upper limit in the integration (in practice it is sufficient to use T up ~ 
(5 — 10 )Trdj as for larger temperatures TP of axinos is more efficient, but the fast expansion of the Universe in 
the reheating period dilutes away all the axinos produced at that early times). We also used A = a/ai = clTr 
and 


d InT 
dA 


J_cLR _ J_ 

4 R d A A 

i , 1 d lng,(T) > 

1 "r" 4 d InT 


(9) 


where R is related to the energy density of radiation pR by R = p_ro 4 . 


The scattering term. In order to deal with the scattering contribution to T~ rp we express S scat in terms of 
the scattering cross-section cr(s), following [SB]), and obtain: 


T/-scat ,i,j 

r a,0 


9i9j 


/■Tup poo 

( d In T \ 

-1 A 2 ] 

/ dT / da- 

/T 0 J (mi+m 2 )/T 

[(’ T cU ) 



soAg 16?r 4 J To 

T 2 Ki(x) <j(x 2 T 2 ) ( x 2 T 2 — m\ — m^Y — 


( 10 ) 


We then change the order of integration and decompose the above integral into 


-w-scat ,i,j 

Y d,0 


1 9i9j 
sqAq 167t 4 



( 11 ) 


where I\ and I 2 are given by expressions very similar to (10), but with T integrated from To to (mi 
and from (mi +rri 2 )/x to T up , respectively. Then / 2 ~ 0 since T 0 


m 2 )/x 

0. For the remaining integral I\ one obtains 


TAScat,i,j 

r a,0 


9 9i 9j Afpi 
167T 4 


rtT,„ 


I (m 1 +m 2 )/T u] 


dt t 3 Ki (t) 


lm 1 +m 2 


d(A0 


f(Vs)cr(s) 


(s — m 2 — m|) 2 — 4 m\m 2 


( 12 ) 
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where g = 135 V 7E and 

y 2 tt3 gl /2 


f(Vs) = 


I 9* 


T 0 3 ^ 3 V 30 


^ -T 


d In Tn - 1 


a 2 t 6 


dA J 


«j> T 4 RT l 

Z_ 4 R I 444 R 

A 3 ~r A 4 


with T = 

t 


(13) 


A careful analysis of Eq. (13) in the reheating period shows that one can make an approximation 


/ = 



(< 1 ) in the reheating period, 
in the RD epoch, 


(14) 


where c ~ — 7 and Trd ~ 0.5 Tr. In practice we find more exact values of a and Trd numerically, but they 
depend only slightly on the model parameters. 


High Tr limit of the scattering term. The integral (12) can be rewritten as a sum of three integrals, 
schematically represented by 


T/'ScatjZ,.? 

r a,0 


Mm 1 +m 2 )/T R r, 

' (m 1 +m 2 )/T up ■ 


rtT m 




rtT m 


+ 


— J\ + + Js- 


m 1 +m 2 J(mi+m2)/%> Jm 1 +m 2 J (m 1 +m 2 )/T r d J tT RD 


(15) 


One can verify that in the limit Trd —> oo, we have J\ —> 0, since the range of external integration shrinks to 
zero , while the integrand does not diverge. The second integral, J 2 , corresponds to the the standard result 
obtained in the instantaneous reheating approximation. In the limit of high reheating temperature the inner 
integral can be simplified to 


rtT-RX) 


d(v / s)/(s)cr(s,t) 


(s — m\ — m ^) 2 — 4 m\m\ 


/ mi+m2 


ptTn d 

~ / d(-\/s) 1 x cr(t) x 1 

J mi+m 2 
~ t a(t) Trd, 


(16) 


(a depends on t via m e g, see, e.g., [333) where we noticed that the integral is mainly determined by the values 
of the integrand in high-s limit in which, to a good approximation, a(s,t ) = cr(t). For the third integral, J 3 , we 
similarly note that inner integration leads to 


f‘ T ~ avD m «(». n *„ (i) r~ d( vs, /( »). 

'tTji d ^'Trd 


(17) 


In the integration range T = y/s/t > T RD and therefore (fl7| becomes 


tcr(t) / dT 
■'Trd 


(W- 


— ta(t) Trd, 
0 


(18) 


where we assumed T up = cTrd with c high enough so that effectively T up can be replaced by 00 in the integration. 
The remaining (external) integrals for both J 2 and J 3 are the same. Hence for high Tr 


J 3 

J '2 


~ 0.17 


Eq. (19) is valid for each contribution to the scattering term. 


(19) 


13 













The decay term. In the case of the decay term we substitute (T) n/ q (see, e.g., [58]) into ^ and obtain 


Y- 


dec.i 1 1 9i TTlj 


a,0 


soblp 27r 2 J To 


dT 

/•OO 

/ dx 

\(-T d \ T ) 

-1 A 2~ 

Too 

J mi/T 

V dA / 

H 


rjnZ 


T 2 


e x =F 1 


( 20 ) 


Once again we change the order of integration and find that one term is negligible, while the other leads to 

r tT u 


v dec,i ggiTnuMpi 

a,° - 2 7r 2 


di ■ 




' m/T u] 


e* =F 1 


d Ef(E)±\l 


E 2 ' 


where / is given by Eq. (14) with s/s replaced by E. The inner integral 




9mi,T R (t) = I d Ef{T = E/t)— i \jl--^, 


to . 


E 2 


( 21 ) 


( 22 ) 


where mi/T up < t < oo can be calculated analytically. Depending on the value of t one obtains for mi/T up < 
t < mi/T RB 


9 mi ,T R (t)=9^ TR (t) = %fC 


6 


3*" 3 - 5” 5 + 7*”’ - 9”” + II" 11 


yi [ m ,./(7V'n 1( )] 


and for t > mi/T RRt (the temperature can be either larger or smaller than Trd) 

9rn t ,T R (t) = g™T R (t) + S^Tr W, 

where ( t R = mi/T RD ) 


9™,T R (t) = 8^3 (I - arctan yp=| + t f^ 2 - ’ 

yi~[mi/(tT up )] 


T 7 f 7 / 1 

4 RD 7 71 


Ofc(*)= m io 


6 7 d 9 ( 1 n 


_ w 3 _ _ w5 + _ w 7 _ _ w 9 + _ w 


sjl- [m;/(tT RD )] 

One can verify that in the case of instantaneous reheating the standard result [55] is rederived. 


(23) 


(24) 


(25) 


(26) 


Appendix C: Phase space integrals for non-thermally produced axinos 


In this appendix we provide results for both the LOSP lifetime and the present-day rms velocity of axinos in 
the case of C 0 yy = 0. Depending on the nature of the LOSP this requires analysis of 3— of 4—body decays. 

3—body bino decay to the axino with C a yY = 0 We calculate the lifetime that corresponds to the 
3—body decay shown in the left panel of Fig. [9] In the case of rria <C mg one obtains 


r ~ 30 sec 


( 100 GeV ) 

5 ( ™S \4 

( 1 TeV ) 

2 

f fa > 

2 

\ m B ) 

V 1 Te V ) 

“s J 

(lO 11 GeV j 

to 

~§ + 3y + (3y 2 - 

T - 1 

+ 

I 1 -!] 

} 


where y = m? /m ~. 


(27) 


It is straightforward to verify that Eq. (27) can be simplified to Eq. (]3 


for to- <C to?. 
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Figure 9: Feynman diagrams for the bino LOSP (left panel) and the stau LOSP (right panel) decaying into the 
axino for C u yy = 0. 


When treating the WDM constraint we calculate the present-day rms velocity 


by 


( v °) = 4.57 x 10 


-5 km 1/12 (|ps|) f T S\ 1/2 


9 d 




(28) 


where (|pa)| is the average momentum of the outgoing axino. The final result reads 

9 


v°d) = 2.5 x 10 


— 3 km\ — 1 / 1.2 ( lGeV\ 


where 


f(y) = 


s ) 

(bs|> 


9 d 


mn 


" b B 

100 GeV 


1 TeV 

m.n 


fa 

10 11 GeV 


x/ 


-1 + 3y + (3y 2 - % + 1) In 


1-i 

v 


~ 1.22 y — 0.47. 


m l\ 

m s) ’ 

(29) 


(30) 


The numerical approximation in the last equation works well unless mg ~ nriq for which function / becomes 
suppressed. 

4—body stau decay to the axino with C & yy = 0 The respective Feynman diagram is shown in the 
right panel of Fig. [9] The final result for the lifetime of the stau reads 


_ , nnri . /100 GeV \ 10 / mg G / m„- \4/lTeV\ 2 / f a \ 

t - l/r{ • sec) x ^ m _ J (lOOGev) (lTev) ( m g ) V10 n GeVy 

rrij 
rria 


X/ 


9 


rrif 

m B 


(31) 


where Cr = 1 for the ’’right” stau and Cl — 4.45 for the ’’left” stau, while functions f(x), g(x) and h(x) are 
shown in Fig. 10 They all tend to 1 for id, i.e., for ma <C m? mg,ra^. 


The presenGday rms velocity is equal to 


km\ _i/i2 


- ( 2.58 x 10“ J 9~d Cl/r x / 


(rnd\„(rnf_\ ~ h (m? 

\ m B 


lGeV^j ^ 100GeV^ 4 ^ 


m~ n 


m-. 


100 GeV 


y ( V 

) VI TeV/ 


2 /1 TeV 


m s 


fa 


10 11 GeV 


(32) 
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Figure 10: Auxiliary functions in the formula for the lifetime of stau decaying into axino if C q yy = 0. 



Figure 11: Auxiliary functions in the formula for the present-day rms velocity of axinos produced in late-time 
stau decays if C a YY = 0. 

where Cr = 1 for the ’’right” stau and Cl — 2.3 for the ’’left” stau while functions /, g and h are shown in 

Fig- [TT] 


Appendix D: Description of the numerical analysis 

Here we explain some details of our numerical analysis of the scenario of axino DM with low reheating tem¬ 
peratures of the Universe in the context of the MSSM. A study of a completely general MSSM would be at 
the same time complicated and unnecessary, hence we select a 10-parameter version of the MSSM (plOMSSM) 
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Parameter 

Range 

bino mass 

0.1 < Mi < 5 

wino mass 

0.1 <M 2 < 6 

gluino mass 

0.7 < M 3 < 10 

stop trilinear coupling 

-12 < A t <12 

stau trilinear coupling 

-12 < A r < 12 

sbottom trilinear coupling 

A b = -0.5 

pseudoscalar mass 

0.2 < rriA < 10 

/i parameter 

0.1 < ji < 6 

3rd gen. soft squark mass 

0.1 < my, < 15 

3rd gen. soft slepton mass 

0.1 < my < 15 

lst/2nd gen. soft squark mass 

top; = Top + 1 TeV 

Wl, 2 W3 

lst/2nd gen. soft slepton mass 

my = my + 100 GeV 

ratio of Higgs doublet VEVs 

2 < tan p < 62 

Nuisance parameter 

Central value, error 

Bottom mass m b (m b ) MS { GeV) 

(4.18, 0.03) [33 

Top pole mass m t (GeV) 

(173.5, 1.0) [51: 


Table 1: The parameters of the plOMSSM and their ranges used in our scan. All masses and trilinear couplings 
are given in TeV, unless indicated otherwise. All the parameters of the model are given at the SUSY breaking 
scale. 


Measurement 

Mean 

Error: exp., theor. 

Ref. 

m h 

125.7 GeV 

0.4 GeV, 3 GeV 


n x h 2 

0.1199 

0.0027, 10% 

m 

BR (B -> X s7 ) x 10 4 

3.43 

0.22, 0.21 

H2 

BR(B U -> r^)xl0 4 

0.72 

0.27, 0.38 

m 

A M Bs 

17.719 ps _1 

0.043 ps" 1 , 2.400 ps" 1 

m 

sin 2 6 cS 

0.23116 

0.00013, 0.00015 

m 

Mw 

80.385 GeV 

0.015 GeV, 0.015 GeV 

ED 

BR (B s ->• n + ir)x 10 9 

2.9 

0.7, 10% 

[Ml Si 


Table 2: The constraints imposed on the parameter spaces of the plOMSSM and the CMSSM. The LUX upper 
limits [70l have been implemented as a hard cut. 


which has practically all the relevant features of the general model. These adjustable parameters of the model 
and their ranges are specified in Table [l] Our choice is closely related to that of m (see discussion therein), 
except that we keep both the wino mass M 2 and the bino mass Mi free. 

We scan the parameter space of plOMSSM following the Bayesian approach. The numerical analysis was 
performed using the BayesFITS package which utilizes Multinest m for sampling the parameter space of 
the model. Mass spectra were calculated with SOFTSUSY-3.4.0 [62], while i?-physics related quantities with 
Superlso v3.3 (Til . 

The constraints imposed in scans are listed in Table [2] The LHC limits for supersymmetric particle masses 
were implemented following the methodology described in Em in]. The DM relic density for low Tr was 
calculated by solving numerically the set of Boltzmann equations, as outlined in [23 HQ. In order to find 
the point where WIMPs freeze out, we adapted the method described, e.g., in [72i to the scenario with a 
low reheating temperature, extracting the relevant functions entering the Boltzmann equations ((<Ti>) e ff and 
(CTu)eff(-E)eff, cf. Ref. m) with appropriately modified MicrOMEGAs v3.6.7 [75] - 
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